#Calculamos la densidad media de conexiones 
# =1 si densidad isotropa
# >1 mas probabilidad que la esperada
# <1 menos probabilidad que la esperada

import math
from pylab import *
import numpy as np


#lista de amplitudes a analizar
cuales=['0.05','0.1','0.2','0.5','1','2','5','10','20','40','60','100','200']
#Directorio de trabajo
dir='/home/jcfernandez/resultados/eps1/u20/AMPL/W2.0/CONNECT/NEO/'
dir='/Users/jonk/Documents/Trabajo/resultados/eps1/u0/AMPL/W2.0/CONNECT/NEO/'
for j in range(len(cuales)):
    #Nombre del fichero a importar
    name='densidadconnectA'+cuales[j]+'.dat'
    name='nconnectA'+cuales[j]+'.dat'
    print 'Fichero: ',name
    #Arrays para el calculo de la media y el error
    media=np.zeros(100)
    media2=np.zeros(100)
    count=np.zeros(100)
    #Importamos los datos
    for line in file(dir+name):
        line=line.replace('\n','')
        val=line.split(' ')
        res=[]
        n=0
        for v in val:
            if v!='':
                res.append(float(v))
        t=res[0]
        densi=res[1::]
        x=[]
        y=[]
        #Discretizacion angular
        step=180./(len(densi))
        for i in range(len(densi)):
            xval=float(i+0.5)*step
            media[i]+=densi[i]
            media2[i]+=densi[i]*densi[i]
            count[i]+=1.0
    
    x=[]
    y=[]
    erry=[]
    #Calculamos los valores  medios y el error
    print len(densi),len(count)
    for i in range(len(media)):
        x.append((i)*step+0.5*step)
        if count[i]>0.0:
            fracArea=(math.cos(math.pi/180.0*(x[i]-0.5*step))-(math.cos(math.pi/180.0*(x[i]+0.5*step))))*0.5
            factor=fracArea*(1000.0*999.0)
            valmedio=media[i]/count[i]/factor
            error=math.sqrt(abs(media2[i]/count[i]-valmedio*valmedio))/factor
            erry.append(error)
            y.append(valmedio)
        else:
            y.append(0.0)
    
    f = open(dir+'dis-'+name,'w')
    for i in range(len(x)):
        f.write(str(x[i])+' '+str(y[i])+' '+str(erry[i])+'\n')
#        f.write(str(x[i])+' '+str(y[i])+'\n')
    f.close()

#plot(x,y,'ro')
errorbar(x, y, yerr=erry, fmt='o-')
show()
